home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / DELPHI32 / MATH / PI / ARCTAN.PAS next >
Pascal/Delphi Source File  |  1996-06-17  |  7KB  |  277 lines
  1. unit arctan;
  2.  
  3. interface
  4.  
  5. uses
  6.   Forms, SysUtils;
  7.  
  8. const
  9.   // If you want you can set these back to 9, 10 and 1 to see the performance
  10.   // with only one digit per array entry.
  11.   MaxValue = 99;
  12.   LessValue = 100;
  13.   MinPlaces = 2;
  14.  
  15. type
  16.   // My initial design called for a base class TArbitraryPrecision that just
  17.   // took parameters (data, value, size). But I decided it would be faster if
  18.   // everything was in the one class as there would be less parameter passing
  19.   // required then. Of course thinking about it later the impact of the
  20.   // parameter passing would have been insignificant compared to the times
  21.   // taken to do the calculations. So perhaps I should have gone for the
  22.   // cleaner design?
  23.   TByte = array[1 .. ((MaxInt div sizeof(Byte)) - 1)] of Byte;
  24.   PTByte = ^TByte;
  25.   TCalculatePi = class
  26.     precision: Integer;
  27.     vm, vn: PTByte;
  28.   private
  29.     procedure SetUp(Divisor: Integer);
  30.     procedure Add;
  31.     procedure Subtract;
  32.     procedure Multiply(Multiplier: Integer);
  33.     procedure Divide(Divisor: Integer);
  34.     procedure WriteSource(FileName: string; Divisor, Iterations: Integer);
  35.     procedure ReadSource(FileName: string);
  36.   public
  37.     constructor Create(SetPrecision: Integer);
  38.     destructor Destroy;
  39.     procedure ArcTan(Divisor: Integer);
  40.     procedure Sum;
  41.   end;
  42.  
  43. implementation
  44.  
  45. uses pidata;
  46.  
  47. constructor TCalculatePi.Create(SetPrecision: Integer);
  48. begin
  49.   PiDataModule.start_position:= 1;
  50.   precision:= SetPrecision;
  51.   GetMem(vm, Precision + 1);
  52.   GetMem(vn, Precision + 1);
  53. end;
  54.  
  55. destructor TCalculatePi.Destroy;
  56. begin
  57.   FreeMem(vn, precision + 1);
  58.   FreeMem(vm, precision + 1);
  59. end;
  60.  
  61. procedure TCalculatePi.SetUp(Divisor: Integer);
  62. var
  63.   index: Integer;
  64. begin
  65.   PiDataModule.start_position:= 1;
  66.   // Calculate vn as 1/K
  67.   fillchar(vn^, precision, 0);
  68.   vn^[1]:= 1;
  69.   Divide(Divisor);
  70.   // And copy it back into vm
  71.   Move(vn^, vm^, precision);
  72. end;
  73.  
  74. // Note: The checks in the while loops in the next three functions on the
  75. // index > 0 is probably superflous
  76. procedure TCalculatePi.Add;
  77. var
  78.   index, carry: Integer;
  79. begin
  80.   carry:= 0;
  81.   index:= precision;
  82.   while (index > 0) and
  83.         ((index >= PiDataModule.start_position) or (carry <> 0)) do
  84.   begin
  85.     vm^[index]:= vm^[index] + vn^[index] + carry;
  86.     if vm^[index] > MaxValue then
  87.     begin
  88.       carry:= 1;
  89.       vm^[index]:= vm^[index] - LessValue;
  90.     end
  91.     else
  92.       carry:= 0;
  93.     Dec(index);
  94.   end;
  95. end;
  96.  
  97. procedure TCalculatePi.Subtract;
  98. var
  99.   index, carry, temp: Integer;
  100. begin
  101.   carry:= 0;
  102.   index:= precision;
  103.   while (index > 0) and
  104.         ((index >= PiDataModule.start_position) or (carry <> 0)) do
  105.   begin
  106.     temp:= vm^[index] - vn^[index] - carry;
  107.     if temp < 0 then
  108.     begin
  109.       carry:= 1;
  110.       vm^[index]:= temp + LessValue;
  111.     end
  112.     else
  113.     begin
  114.       carry:= 0;
  115.       vm^[index]:= temp;
  116.     end;
  117.     Dec(index);
  118.   end;
  119. end;
  120.  
  121. procedure TCalculatePi.Multiply(Multiplier: Integer);
  122. var
  123.   index, carry, temp: Integer;
  124. begin
  125.   carry:= 0;
  126.   index:= precision;
  127.   while (index > 0) and
  128.         ((index >= PiDataModule.start_position) or (carry <> 0)) do
  129.   begin
  130.     temp:= vn^[index] * Multiplier + carry;
  131.     carry:= temp div LessValue;
  132.     vn^[index]:= temp mod LessValue;
  133.     Dec(index);
  134.   end;
  135. end;
  136.  
  137. procedure TCalculatePi.Divide(Divisor: Integer);
  138. var
  139.   index, carry, temp: Integer;
  140. begin
  141.   carry:= 0;
  142.   for index:= PiDataModule.start_position to precision do
  143.   begin
  144.     temp:= vn^[index] + (carry * LessValue);
  145.     carry:= temp mod Divisor;
  146.     vn^[index]:= temp div Divisor;
  147.   end;
  148. end;
  149.  
  150. procedure TCalculatePi.ArcTan(Divisor: Integer);
  151. var
  152.   iterations, square, check_zero, index: Integer;
  153.   do_add: Boolean;
  154. begin
  155.   // Initialise all the local variables
  156.   iterations:= 0;
  157.   do_add:= True;
  158.   square:= Divisor * Divisor;
  159.   // Initialise the vn'th array element to 1/Divisor and copy it to the vm'th
  160.   // element
  161.   SetUp(Divisor);
  162.   // Now iterate until we have zeroed the vn array
  163.   // Note: That there is always an error in the last few digits due to not
  164.   // having enough precision. Hence to get n digits you might like to do
  165.   // n+50 or some such
  166.   while PiDataModule.start_position < precision do
  167.   begin
  168.     do_add:= not do_add;
  169.     // Note the following divide needs to be done in two steps as the value
  170.     // ((2 * iterations) + 3) * square can exceed MaxInt on large runs
  171.     Divide(2 * iterations + 3);
  172.     Divide(square);
  173.     Multiply(2 * iterations + 1);
  174.     if do_add then
  175.       Add
  176.     else
  177.       Subtract;
  178.     // Calculate the new location of the first zero in the vn array
  179.     check_zero:= PiDataModule.start_position;
  180.     while (vn^[check_zero + 1] = 0) and (check_zero < precision) do
  181.     begin
  182.       Inc(PiDataModule.start_position);
  183.       Inc(check_zero);
  184.     end;
  185.     Inc(Iterations);
  186.     Application.ProcessMessages;
  187.   end;
  188.   // And now write the data to a storage file
  189.   WriteSource('output' + IntToStr(Divisor) + '.txt', Divisor, iterations);
  190. end;
  191.  
  192. // Yes, I admit it, I was incredibly lazy here, but it works.
  193. procedure TCalculatePi.WriteSource(FileName: string;
  194.                                    Divisor, Iterations: Integer);
  195. var
  196.   index: Integer;
  197.   temp: string;
  198.   output_file: Text;
  199. begin
  200.   assign(output_file, FileName);
  201.   try
  202.     rewrite(output_file);
  203.     if (Divisor = 0) and (Iterations = 0) then
  204.     begin
  205.       writeln(output_file, 'PI = 3.+');
  206.       writeln(output_file);
  207.       for index:= 2 to precision do
  208.       begin
  209.         temp:= IntToStr(vm^[index]);
  210.         while Length(temp) < MinPlaces do
  211.           temp:= '0' + temp;
  212.         write(output_file, temp);
  213.         if (index - 1) mod 100 = 0 then
  214.           writeln(output_file)
  215.         else if (index - 1) mod 10 = 0 then
  216.           write(output_file, ' ');
  217.       end;
  218.     end
  219.     else
  220.     begin
  221.       writeln(output_file, 'Arctan(1/' + IntToStr(Divisor) + ')');
  222.       writeln(output_file);
  223.       writeln(output_file, 'Iterations = ' + IntToStr(Iterations));
  224.       writeln(output_file);
  225.       for index:= 1 to precision do
  226.         writeln(output_file, IntToStr(vm^[index]));
  227.     end;
  228.   finally
  229.     close(output_file);
  230.   end;
  231. end;
  232.  
  233. procedure TCalculatePi.ReadSource(FileName: string);
  234. var
  235.   index: Integer;
  236.   input_file: Text;
  237.   a_bit: string;
  238. begin
  239.   assign(input_file, FileName);
  240.   try
  241.     reset(input_file);
  242.     readln(input_file);
  243.     readln(input_file);
  244.     readln(input_file);
  245.     readln(input_file);
  246.     for index:= 1 to precision do
  247.     begin
  248.       readln(input_file, a_bit);
  249.       vn^[index]:= StrToInt(a_bit);
  250.     end;
  251.   finally
  252.     close(input_file);
  253.   end;
  254. end;
  255.  
  256. procedure TCalculatePi.Sum;
  257. var
  258.   output_file: Text;
  259. begin
  260.   PiDataModule.start_position:= 1;
  261.   ReadSource('output18.txt');
  262.   Multiply(48);
  263.   Move(vn^, vm^, precision);
  264.   PiDataModule.start_position:= 1;
  265.   ReadSource('output57.txt');
  266.   Multiply(32);
  267.   Add;
  268.   PiDataModule.start_position:= 1;
  269.   ReadSource('output239.txt');
  270.   Multiply(20);
  271.   Subtract;
  272.   WriteSource('pi.txt', 0, 0);
  273. end;
  274.  
  275. end.
  276.  
  277.